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We consider two fully frustrated Ising models: the antiferromagnetic triangular model in a field 
of strengtli, h = HTkn, as well as the Villain model on the square lattice. After a quench from a 
disordered initial state to T = we study the nonequilibrium dynamics of both models by Monte 
Carlo simulations. In a finite system of linear size, L, we define and measure sample dependent 
"first passage time", tr, which is the number of Monte Carlo steps until the energy is relaxed to the 
ground-state value. The distribution of tr, in particular its mean value, {tr{L)), is shown to obey the 
scaling relation, {tr{L)) ~ \n{L/Lo), for both models. Scaling of the autocorrelation function of 
the antiferromagnetic triangular model is shown to involve logarithmic corrections, both at H = 
and at the field-induced Kosterlitz-Thouless transition, however the autocorrelation exponent is 
found to be H dependent. 



I. INTRODUCTION 

The phenomenon of ageing in physics is first noticed in 
glassy systemsi— , when the physical properties of materi- 
als are found to depend on the preparation process. In a 
microscopic point of view during relaxation these systems 
goes over an infinite set of mctastablc states and thermo- 
dynamic equilibrium is never reached in finite time. Later 
it was observed that similar phenomenon could take place 
in simple homogeneous ferromagncta^, too, if the the sys- 
tem is quenched from the paramagnetic phase below or at 
the critical temperature, Tc. After the quench domains 
of parallel spins are started to grow and the typical size 
of a domain of correlated sites, £,{t), after time t is given 

by 



(1) 



z being the dynamical exponent. In the case of a critical 
quench generally z corresponds to the dynamical expo- 
nent measured at equilibrium and it generally does not 
depend on the temperature of the initial state, Ti. This 
is, however not true for the two-dimensional XY model 
in which case at the two sides of the Kosterlitz-Thouless 
temperature, Tkt, we have^: 



(t/lnOi/2, 



if Ti < Tkt, 
if T, > Tkt. 



(2) 



Due to the quench the translational symmetry in time 
is broken consequently the autocorrelation function, 
G{t,.s), measured after a waiting time s < t is non- 
stationary. In the limit s ^ 1 it is expected to behave 



G(i,s) 



(3) 



whereas for the XY model the scaling function depends 
on the combination^, y — tins/ slnt: 



G{t,s) 
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In both cases the exponent, a, is related to the decay 
exponent of critical spatial correlations: 

d-2 



V 
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and the scaling function, G{y), for large arguments be- 
haves as^: 



G{y) y-V^ 



(6) 



where the autocorrelation exponent, A, is a new param- 
eter which can not be expressed in terms of other (equi- 
librium) exponents. 

There is a class of systems in which due to competing 
interactions or frustration the critical point is shifted to 
Tc = 0, thus there is no ordered phase. The prototype of 
such systems are the antiferromagnetic triangular Ising 
modcl^ii^ (ATIM) and the Villain modelii, which is the 
fully frustrated Ising model (FFIM) on the square lat- 
tice. Both models have a finite ground state entropy per 
site and the usual ensemble average for T > is replaced 
by averaging over the different degenerate ground state 
configurations. It is known exactly for several decades 
that the ATIM and the FFIM are in the same univer- 
sality classi^'i^'i^ii^ as far as the decay of the spin-spin 
correlation function is considered. In the presence of a 
weak magnetic field of strength, h — HTks, the ATIM 
is shown to stay criticali^ for H < He. At the end of 
the critical line at H = He there is a Kosterlitz-Thouless 
transitionii, which is analogous to that observed in the 
XY model. 

As far as the nonequilibrium critical properties of the 
ATIM and the FFIM after a quench to Tc = are con- 
cerned they are also expected to be in the same univer- 
sality class, but about the value of the exponents, such 
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as z and A, there are conflicting results in the literature, 
which all refer to -ff = 0. In RefJ^ the authors suggest 
that in the ATIM there are topological defects which are 
attracted or repelled by entropy-driven Coulomb forces, 
which are analogous to vortices and anti-vortices in the 
two-dimensional XY model for T > Tkt- Therefore 
noncquilibrium relaxation in the two models arc expected 
to be governed by the same mechanism and noncquilib- 
rium scaling is given in the same form, see the second 
equation of Indeed numerical calculation of the den- 
sity of topological defects are in favour of these consid- 
erations. On the other hand, in Refi^^ the numerical 
data are found in favour of a distinct dynamical expo- 
nent, z ~ 2.33, which is explained in terms of interaction 
between defects and loose spins. In this study the au- 
tocorrelation function is found to involve the exponent, 
\/z ~ 0.86. More recently the noncquilibrium dynam- 
ics of the FFIM after a quench to Tc = is studied 
numericall}*2S and the obtained data are interpreted in 
favour of the same type of scaling as in the XY model 
for T > Tj^T, with z = 2 supplemented by logarith- 
mic corrections. 

We note on studies of related models which can be 
obtained through mappings to two-dimensional cover- 
ing problems^ii^, such as by rhombi or dominoes or by 
dimers in the dual lattices. As described in section |TT] 
these are also equivalent to solid-on-solid (SOS) mod- 
els. Recently it has been shown^^ that the p-terphenyl- 
3,5,3',5'-tetracarboxylic acid adsorbed on graphite self- 
assembles into a two-dimensional rhombus random tiling. 
Its relaxation dynamics and relation with the similar pro- 
cess in structural glasses is studied in Re^. The relax- 
ation dynamics of the compressible ATIM, in which also 
elastic strain fields are taken into account is studied in 
Refi^, whereas the ATIM in a staggered magnetic field 
is considered in Re^^. 

In this paper wc study the noncquilibrium critical dy- 
namics of the ATIM and the FFIM by Monte Carlo 
(MC) simulations and try to solve the existing contro- 
versy about their scaling properties. Our approach dif- 
fers from the previous ones in two aspects. First, in a 
finite system of linear size, L, we define and measure 
a quantity, the sample (or initial condition) dependent 
first passage time, tr. Its distribution and in particular 
its mean value, {tr{L)), is used to obtain a direct relation 
between the characteristic time and the length in the sys- 
tem and through Eq.(IT]) we obtain a direct estimate for 
the dynamical exponent. We note that in previous nu- 
merical studies the calculated critical quantities are char- 
acterised by combination of different critical exponents, 
see e.g. the autocorrelation function in Eq.Q, and no di- 
rect estimate for z has been made. The second feature of 
our study is that for the ATIM we consider not only the 
zero field model, but the Kosterlitz-Thouless transition 
at H = He as well. We compare the numerical results in 
the two cases in order to decide if the two transitions are 
in the same dynamical universality class. 

The structure of the paper is the following. In SccHIl 



the models are described. In Sec lIIII wc present results 
about the sample dependent first passage times both for 
the ATIM (at i? = and at H ^ He) and for the FFIM. 
The autocorrelation function of the ATIM is studied in 
SeclIVI We close our paper by a discussion. 

II. MODELS 

We consider two fully frustrated Ising models. The 
first model is the antiferromagnetic Ising model on the 
triangular lattice in the presence of an external field, h, 
which is in the order of the thermal fluctuations, fc^T. 
The Hamiltonian is defined in terms of = ±1 as 

^^-/vJ^a.a.-H^a,, (7) 

where i and j are sites of a triangular lattice and the first 
sum runs over nearest neighbours. We consider the limit 
T ^ 0, so that K — > — oo but at the same time wc have 
H = 0{1). 

The model is exactly solved at zero field, when the 
entropy per site is^iiiS: 

S^ - [ ' ln(2 cos a;)dw ~ 0.32306 , (8) 

which is calculated for periodic boundary conditions. We 
note that the actual form of the boundary conditions is 
crucial for the ground state entropy. For example con- 
sidering a hexagonal part of the lattice with open bound- 
ary conditions the entropy can be any value in the range 
[O, I ln3 — 21n2] depending on the ratio of the lengths 
of the sample^!. On the other hand for a part with paral- 
lelogram shape and open boundary conditions there are 
only two ground states and consequently the entropy is 
zero2^. The spin-spin correlation function at = is 
calculated by Stephcnsor>i^ which shows algebraic decay: 
thus 77 = 1/2. 

Further interesting results are obtained through dif- 
ferent mappings, which are valid also for non-equal cou- 
plings as well as in the presence of external fields. The 
model is equivalent to a covering problem with rhombi in 
the triangular lattice or with dimer covering in the dual 
hexagonal lattice^ii^. Assigning appropriate height vari- 
ables to the rhombi one obtains a solid-on-solid (SOS) 
model which has a roughening transition at H ^ He 
of the Kosterlitz-Thouless type and the rough phase is 
present for small fields. The SOS model is then trans- 
formed by a standard although non-exact mapping to a 
Gaussian mode l^^'^^ from which the conformal anomaly 
is obtained as c = 1 and different spin- and vertex- 
exponents are calculated at H = 0. For non-zero 
field the model is studied numerically by transfer-matrix 
techniquesii. The system is shown to be critical for 
< H < He and has long range order otherwise. The 
location of the Kosterlitz-Thouless transition is found at 
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He — 0.266. Some exponents arc found to be constant 
along the critical line, whereas others are 7?-dependent 
but they satisfy the weak-universality hypothesis, i.e. the 
ratio of specific exponents is constant. At the Kosterlitz- 
Thouless transition point the spin-spin correlation func- 
tion is conjectured to decay with an exponent: rjc = 4/9. 

The second model is the fully frustrated Ising model 
on the square lattice defined by the Hamiltonian: 

^ = -I^^^^.j^.f^j- (9) 

In the even vertical lines the nearest neighbour couplings 
are antiferromagnetic, Kij ~ ^K, otherwise they are 
ferromagnetic, Ki j = K . Consequently at each elemen- 
tary square the product of the couplings have a negative 
sign, and the external field is = 0. This model can 
be mapped to a domino covering problem in the square 
latticoii or equivalently to a dimcr covering problem in 
the dual lattic o^^i'^° . From this mapping one can obtain 
the ground state entropy per sit o^^i^^'^^ : 

5 = - = 0.2916 (10) 

TT 

where G is the Catalan's constant. The spin-spin cor- 
relation function is calculated cither through a map- 
ping to the eight-vertex modcU^ or by transfer matrix 
techniqucsi^ii^. The decay exponent is found to be 
rj ~ 1/2, which is the same as for the ATIM at _ff = 0. 
In a similar way as the ATIM, the FFIM is mapped to a 
SOS model on the square lattice^. Analysing the relax- 
ation process of this model in equilibrium a lower bound 
for the relaxation time is obtained: t{L) > 0{L^), thus 
z > 2. This relation in the continuum approximation 
turns to an equality. 

III. RELAXATION TIME IN FINITE SYSTEMS 

We study the models by MC simulations. They are 
prepared in a random initial state (i.e. with T = oo) 
and then quenched to Tc = and let to evolve ac- 
cording to the Glauber dynamics. A randomly cho- 
sen spin at site i is flipped with a probability: p{i) = 
^-AEUksT i^^AEjksT ^ g-AB./fc^T)^ ^^ere the reduced 

energy-difference is: —AEi/ksT = (Ti{K Y^'- cTj+H) and 
the sum in j runs over nearest neighbours of i. In the 
limit we consider (A^ —oo) a flip is always accepted if 
the interaction energy is lowered (1** type flip) and ac- 
cepted with a probability: pi = e~^°'' /{e^°'' + e~^'^^), 
for loose spins, i.e. for which aj ~ (2"'' type flip). 

During the relaxation process in an infinite system the 
interaction energy is monotonously decreasing due to P* 
type flips. In a finite sample, however, at time tr the 
interaction energy reaches its minimum, corresponding 
to a ground-state configuration of the H — problem. 
For later steps, t > tr, only 2"'^ type flips are possible 
and the system evolves in the ground-state sector, i.e. 



the dynamics is in equilibrium. Consequently for a given 
initial state (sample) tr is a characteristic time that we 
call a first passage time, and we assume that it scales 
as a usual relaxation time. Here we use periodic bound- 
ary conditions. For open boundaries possibly one should 
modify the definition of i,., c.f. measuring the time until 
the energy is lowered by a given fraction, but the same 
qualitative picture is expected as for periodic lattices. 

In the following of this section we study the distribu- 
tion and properties of tr for different finite sizes, L. Note 
that the Monte-Carlo algorithm is not ergodic and does 
not sample properly the ground state properties'^, it is 
merely a single spin flip dynamic. We expect that the 
results we present below will remain unchanged if we use 
another dynamic, as soon as the number of spins which 
can be flip at a time is kept fixed when L varies. 

In the actual calculation, as already mentioned in sec- 
tionini we consider finite systems with sites, with peri- 
odic boundary conditions and with L, which is a multiple 
of 3. We have checked that the asymptotic dependence 
of tr with L is the same if L = 3/ -I- 1 or L = 3^ -I- 2, but 
with an increased prefactor and therefore the computa- 
tion time is also larger. The ground states arc usually 
classified using the notion of strings (see Ref i^) , but the 
dynamic always drives the system into the most probable 
string sector. 

We went up to L = 768 and the number of realiza- 
tions are given in Table |T] for the different models. As an 
example the accumulated distribution of tr, Phitr), for 
L = 192 is shown in Fig[T]for the ATIM at = 0. 




FIG. 1: (Color online) Accumulated distribution of the sam- 
ple dependent first passage time, tr, for the ATIM at _ff = 
with linear size, L = 192. For large arguments the tail of the 
distribution is exponential, as shown in the inset. 

The tail of PL{tr) for large argument is exponential, 
Phitr) ~ exp(— ir/T(L)), as shown in the inset of Fig[T] 
in a lin-log plot. For a given size, L, we have calculated 
the mean value, {tr{L)), which is presented in Table [J 
We have checked, that the ratio {tr{L)) /t{L) is approx- 
imately independent of L, thus there is just one time- 
scale in the problem. To see the size-dependence of the 



time-scale we have plotted in Figl2] {tr{L)) vs. L in a 
double-logarithmic scale for the three different problems 
we studied. 
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FIG. 2: (Color online) Average first passage time, {tr{L)), 
as a function of L in a log-log plot for the three frustrated 
models studied in the paper. The straight full line has a slope 
z' = 2.30 obtained by simple linear fit. Note that for large L 
the local slope between the points is smaller than z' , see also 
Table ID 



We have also explicitly confronted the two possible sce- 



Ti(L) = AL" 



(12) 
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FIG. 3: (Color online) Fitting the measured average first pas- 
sage times, (tr (-£/)), of the ATIM at = by the two possible 
functional forms in Eas.(|12p and (|13p . In the inset extrapola- 
tion of the effective dynamical exponents in Table |l] are given 
as a function of l/ln(L/Lo) with Lq = 2.5. 



In all cases the points are approximately on parallel 
straight lines, which would correspond to the relation in 
Eq.([T]) by substituting ^ by L. Having a simple linear 
fit through the points we obtain z' « 2.30, which is in 
agreement with the estimate of z = 2.33 found in Ref 
for the ATIM with H = 0. However, having a closer look 
at the points in Fig[2] on can see that the local slopes 
of the curve is monotonously decreasing with increasing 
L. We have calculated effective exponents, Zeff{L), from 
two-point fit comparing results with L and L': 



effiL) = 



\n{U{L))-ln{U{L')) 
ln(L) - In(i') 



(11) 



which are presented in Table [H too. 

For a given size, L, within the error of the calculation 
the Zeff{L)-s are found to be the same for the three prob- 
lems and they are decreasing below z' — 2.30. Wc have 
extrapolated them by assuming a logarithmic correction 
in the form of 1 / \n{L / Lq) , where Lq is an 0(1) constant. 

This extrapolation is shown in the inset of FigO for 
the ATIM at = 0. The effective exponents for dif- 
ferent L lay on a straight line of slope 1 and the ex- 
trapolated value is close to 2, which is in agreement 
with the scenario in Reffi^, i.e. logarithmic corrections 
to the z = 2 leading behaviour, as given for the XY- 
model in Eq. ^ . Indeed inverting the second equation of 
([2]) we have (if ^ is replaced by L) tr ~ L^ln(L), from 
which the effective exponent is given in leading order as: 
Ze//(i)^2 + l/ln(L/Lo). 



and 



T2(L) = ALMn(L/Lo) 



(13) 



with the data {tr{L)) in Table [J We have minimised the 
quantity 



E 



rfc(L) - (t,(L)) 



rfe(L) + {tr[L)) 



(14) 



for k = \ and 2 and in the sum L runs over the con- 
sidered sizes. For the ATIM with £f = the best fit 
in the two cases are obtained with: ti{L) = O.ISL^'^^, 
and T2{L) = 0.183L^ ln(L/3.6), which arc shown in the 
main part of Fig|31 Indeed, the second scenario yields 
a much better fit and note that in this case the value 
of the dynamical exponent, z = 2, is obtained from the 
optimisation process, its value was not fixed in advance. 
The value of the quantity in Eg. p^ is slowly varying 
if z is in the range [1.96,2.04] with a minimum located 
precisely at z = 2, outside this range, we observe a rapid 
increase. Repeating this analysis for the two other frus- 
trated models we have obtained the same conclusion. For 
the ATIM atH^Hc the best fits are ri(L) = O.ISL^-^^, 
and T2(L) = 0.20L2 ln(L/3.9), whereas for the FFIM 
Ti{L) = 0.20L2-27^ and r2(L) = 0.222.2 ln(L/4.1). 
both cases the relative error of the fits were 5 to 10-times 
smaller for the second scenario. 
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TABLE I: Average of the sample dependent first passage time (tr) for different finite systems of tlie ATIM at zero field and at 
the field induced Kosterlitz-Thouless transition point, H — He, as well as for the FFIM at zero field. The effective dynamical 
exponent calculated from two consecutive sizes, see Eq. pif) is also indicated together with the number of samples used in the 
simulation. 
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IV. AUTOCORRELATION FUNCTION OF THE 
ATIM 

For the FFIM the scahng behaviour of the autocorre- 
lation function in the ageing regime has been thoroughly 
studied in Rcfj2£. In that work the numerical data for 
G{t, s) has been analysed by assuming two possible scal- 
ing forms, as given in Eqs.Q and (|4]) and a clear evidence 
in favour of the second scenario is found. To summarise 
for the FFIM the scaling function is suggested to contain 
logarithmic corrections and the autocorrelation exponent 
is measured X = d. Here we repeat this type of analysis 
for the AFIM, both at 77 = and at iJ Hr. 
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FIG. 4: (Color online) The equilibrium autocorrelation func- 
tion, Geq{t, s) as a function of t — s for the ATIM at = in 
a log- log plot. The asymptotic slope of the curve is given 
by a{H = 0) = 0.2501(10). Inset: the same analysis at 
H = He- The asymptotic slope of the curve is given by 
a(H = He) = 0.229(5). 




FIG. 5: (Color online) Scaling of the reduced autocorrelation 
function, G(t, s) = s^Git, s) of the ATIM at = with 
a = 0.25 assuming conventional scaling, see Eq.(|3}. 



First we consider the zero field model and analyse the 
autocorrelation function at equilibrium, i.e. formally the 
measurement is performed for t,s > tr- In the actual 
calculation we have started from a ground state and after 
thermalization we have performed the measurement. In 
that case due to time translational invariance Geq{t, 5) is 
stationary and we have asymptotically (i.e. for t—s 3> 1): 



Ge,(i,s) w (<-s)- 



(15) 



This behavior is illustrated in Fig. 3] in which the de- 
cay exponent is measured as a{H = 0) = 0.2501(10), in 
agreement with the scaling relation in Eq.([5l). 

Next we consider the autocorrelation function in the 
ageing regime, i.e. for t, s < tr and analyse the scaling 
function, G{t,s) = s°'G{t,s), assuming the two scenar- 
ios in Eqs.Q and (|4]). For the decay exponent we used 
the presumably exact value, a ~ 0.25. For conventional 
scaling, i.e. with y ~ t/s, as given in Eq.Q the results 
are shown in FigO , whereas for ATK-type scaling., i.e. 
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FIG. 6: (Color online) The same as in Fig. [5] with XY-type 
scaling, see Eq.(|4]). 



with y = tins /slut, as given in Eq.Q the results are 
in Figl6l As seen in FiglS] for conventional scaling the 
collapse of the curves for different values of the waiting 
times, s, is certainly not perfect. In the enlarged part 
of the figure presented in the inset is clearly seen a shift 
between the curves and this shift can be approximately 
characterised by Ins. On the other hand in Figl6] for 
XY-type scaling the collapse of the curves is much bet- 
ter: for the three largest waiting times: s = 160, 640 and 
1280 the curves practically fall into each other. Conse- 
quently our data arc in favour of logarithmic corrections 
in the scaling function of the autocorrelation function of 
the ATIM at iJ = 0. We have then estimated the au- 
tocorrelation exponent, A/z, from the slope of the linear 
part of Fig|B] for the three largest waiting times and we 
obtained: X/z = 1.003(10). That value corresponds to 
A = d, as obtained for the FFIM^S. 
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s=1280 



1 2 3 4 5 6 7 

ln(t/s) 

FIG. 7; (Color online) Scaling of the reduced autocorrelation 
function, G{t,s) = s"G(t,s) of the ATIM at H ^ He with 
a = 2/9 assuming conventional scaling, see Eq.Q. 
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FIG. 8: (Color online) The same as in Fig. Qwith XY-type 
scaling, see Eq.Q. 



B. ATIM at H = He 

We have repeated the analysis of the autocorrelation 
function at the field-induced Kosterlitz-Thouless transi- 
tion point, at H = He- In that case the equilibrium 
autocorrelation function is shown in the inset of Fig|4l 
in which case we have measured the decay exponent as 
a{H ~ He) = 0.229(5), close to the conjecturedly exact 
value: rjc/z = 2/9. Results about the autocorrelation 
scaling function are shown in Fig[7]for conventional scal- 
ing and in Fig[5] for XF-type scaling. As for the zero 
field case the scaling collapse of the curves for different 
values of the waiting time are much better for XY-type 
scaling. We have also measured the autocorrelation ex- 
ponent with the result: X/z = 0.93(2). This is somewhat 
smaller than those found in the other two models. 



V. DISCUSSION 

We have studied the nonequilibrium dynamics of dif- 
ferent fully frustrated Ising models in which the critical 
point is at r = 0. Two of these models, the ATIM at zero 
field and the FFIM are in the same universality class, as 
far as equilibrium critical properties are concerned, for 
example the decay exponent is = 1/2 for both models. 
The third model, the ATIM at the critical field, H = 
has a Kosterlitz-Thouless transition and also the correla- 
tion exponent is different. We have quenched the models 
from a disordered initial state to the critical temperature, 
i.e. to T = 0, and studied the nonequilibrium domain 
growth process as well as the scaling properties of the 
autocorrelation function. 

The domain growth is studied at a fixed size of the 
system, L, by measuring the time, tr, until equilibrium 
is reached. Although tr varies from sample to sample 
its distribution is characterised by only a single time- 
scale, which can be taken as the average value , (tr). 
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For all the three models {tr) is found to have the same 
asymptotic L dependence, which is in the form known for 
the XY- model in Eq. ([2|) . Since both the XF-model and 
the ATIM aX H = He have the same Kosterlitz-Thouless 
transition this equivalence found in the non-equilibrium 
critical dynamics is not surprising. 

Using this result we can estimate the asymptotic be- 
haviour of the density of defects, in the systems. To 
be concrete we consider the ATIM at if = for which 
a defect is represented by an elementary triangle having 
all the three spins in the same state. In a finite sample of 
linear size, L, and sample dependent first passage time, 
tr, at t = tr there is one (or 0(1)) defect, consequently 
the density of defects sA, t = {tr{L)) is p{t) ~ L^"^ . Us- 
ing the relation in Eq. ([2|) and replacing ^ by i we arrive 
to the result: p{t) ~ log(t/to)/i, in agreement with the 
conjecture in Ref.— . 

The analogy with the nonequilibrium scaling of the 
XF-model is found to be valid for the autocorrelation 
function, too. For the ATIM, both at zero field and at 
H = He optimal scaling collapse of the autocorrelation 



function is found if the scaling variable is taken in the 
same form, as in the XF-model. For the ATIM at if = 
and in the FFIM the autocorrelation exponent is found 
to be A = d, like in mean-field theory, however with log- 
arithmic corrections. On the other hand in the ATIM 
at the field-induced Kosterlitz-Thouless transition point 
the autocorrelation exponent is found somewhat smaller, 
than d. Consequently in the ATIM with varying field 
the dynamical exponent is presumably constant, z = 2, 
the scaling functions contain same type of logarithmic 
corrections, but the autocorrelation function weakly de- 
pends on the value of H . 
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